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We investigate the flow in a spherical shell subject to a time harmonic oscillation of 
its rotation rate, also called longitudinal libration, when the oscillation frequency is 
larger than twice the mean rotation rate. In this frequency regime, no inertial waves 
are directly excited by harmonic forcing. We show however that it can generate through 
non-linear interactions in the Ekman layers a strong mean zonal flow in the interior. An 
analytical theory is developed using a perturbative approach in the limit of small libration 
amplitude e and small Ekman number E. The mean flow is found to be at leading order 
an azimuthal flow which scales as the square of the libration amplitude and only depends 
on the cylindrical-radius coordinate. The mean flow also exhibits a discontinuity across 
the cylinder tangent to the inner sphere. We show that this discontinuity can be smoothed 
through multi-scale Stewartson layers. The mean flow is also found to possess a weak 
axial flow which scales as 0(e 2 E 5 ^ 42 ) in the Stewartson layers. The analytical solution is 
compared to axisymmetric numerical simulations and a good agreement is demonstrated. 
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1. Introduction 

Most astrophysical bodies, in addition to their rotation, are subject to harmonic forcing 
such as precession, tidal deformation and latitudinal/longitudinal librations. The nonlin- 
ear corrections induced by such a forcing can lead to important steady and axisymmetric 
flow, called zonal flow, in liquid core of telluric planets, in subsurface ocean of some satel- 
lites or in atmosphere of gaseous planets. In this study, we focus on longitudinal libration 
which corresponds to an harmonic oscillation of the rotation rate. These oscillations are 
generally induced by gravitational coupling be tween an astrophysical b ody and its main 
gravitational partner around which it orbits ( Comstock fc Billsl [2003) . The librational 



forcing can have a non-negligible contribution in the internal dynamic of a planet and 
a better knowledge of this dynamic can bring important information on the structure 
of a planet, as for example an estimation of the thickness of the fluid layer under the 
ice shell in some satellites like Europe, Ganymede or Encelade (ISpohn fc Schubertl 2003 



Lorenz et al. 2008 : Van Hoolst et al. 2008; iRambaux et a.Z.1 l201lh or the existence of an 



inner liquid core in telluric planets like Mercury ([Margot et aZ.M2007t ) . 

In the present work, we derive rigorous results for fast longitudinal libration when the 
libration frequency is at least twice larger than the spin frequency. Such a configuration is 
not the most common situation as the main libration frequency often corresponds to the 
orbital frequency which is generally equal or smaller than the spin frequency. Here are a 
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few situations which could exhibit fast libration. For example, one could imagine a planet 
or a satellite for which spin and orbital f requencies are in resonance with a ratio 1:2. This 



situation which is in principle possible ( Goldreich fe Peale 1966 : Wieczorek et al. 2012f) 



would lead to fast libration with a libration frequency equal to twice the orbital frequency. 
Fast libration could also be present in a spin-orbit synchronized object when the forcing 
is induced by another satellite or planet with a larger orbital frequency. The Galilean 
satellites (Io, Europa, Ganymede) are in this situation. Io has an orbital frequency twice 
larger than Europa, and four times larger than Ganymede. It could a priori then force the 
libration of Europa and of Ganymede with a frequency twice and four times larger than 
their spin frequency respectively. In exoplanet systems, a similar situation can also be 
found. For instance, the two first planets in the 55 Cancri system have orbital frequencies 
of 2.81 and 14.65 days ( Fischer et al. 20081 ). Since both planets are expected to be spin- 
orbit synchronized, the frequency of libration of the second planet induced by the motion 
of the first planet would be ui = 14.65/2.81 « 5.2. 

The firs t analytical study of nonlinear effects generated by harmonic forcing was per- 
formed bv lBussd (| 19681 ) in the case of precession in a spheroid. He showed that a zonal flow 
is gene rated in the b ulk as an effect of the nonlinear interactions in the viscous boundary 
layers. IWand (|l970h studied theoretically and experimentally the nonlinear effect s asso- 
ciated with an oscillation of the rotation rate in a cylindrical container (se e also IBussel 
(1201061) 1. The case of a fixed tidal deformation was considered bv lSuessl (Il97ll) who argued 
that a counter-rotating vortex should exist near the axis of rotation. lAldridge fc Toomre 
(1969) analyzed librating spheres but their experimental study focused on the linear re- 
sponse in the presence of inertial waves: they showed using pressure measurements that 
incrtial waves can be resonantly excited in a libr ating sphere for particular libr ation fre- 
quenc ies. This has been confirmed numerically bv lRieutord ( 1991 ). In his thesis, Aldridge 
(|l967l) also mentioned that when the inertial waves are not significant, the flow in the 
bulk seems to drift with an amplitude proportional to the s quare of t he lib ration ampli- 
tude but he did not report any quantitative measurement. iTilgner (1999) investigated 
numerically the linear response to the librational forcing in a spherical shell and studied 
the effect of the presence of an inner core but he focused on the inertial waves and the 
attractors only. Recently, ICalkins et al. (l2010h studied numerically the mean zonal flow 
in the same geometry. They confirmed the presence of a mean zonal flow independent of 
the Ekman number and which scales as the square of the libration amplitude. However, 
their study remained mainly limited to a fixed fr equency. Some LD V measurements of 
the mean zonal flow have also been pe r forme d by iNoir et all (|2010l ) in a cylindrical ge- 
ometry. In this geometry. ISauret et al. ( 20121 ) analyzed numerically the influence of the 
librati on frequency . They found that the mean zonal flow was well described by the the- 
ory of lWangl ( 19701 ) as long as no inertial waves are excited. The on l y theor etical study of 
the zonal flow in a spherical geometry is the recent work bv IBussel ( 2010a ). He obtained 
by asymptotic methods an analytical prediction of the mean zonal flow in the limits of 
small libration frequencies and small Ekman numbers. The present paper extends this 
theory to the case of a spherical shell for any frequency. 

The study implicitly assumes the absence of instability. Both shear and centrifugal 
instability are a priori possible when the amplitude of libration becomes sufficiently 
large. Centrifug a l inst abilit y has been analyse d in details in a cylindrical geometry by 
lErn fc Wesfreidl (Il999h and I Sauret et all (|2012h . It has also been observed in a librating 
sphere by Noir et all (|2009[ ). Shear instability is present in planar Stokes layers (e.g. 
Davisl 1976| ) but no clear evidence has been provided so far in cylindrical and spherical 
geometries. In the present study, the forcing amplitude is assumed sufficiently small 
such that the flow generated at the libration frequency remains stable. To fully resolve 
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this leading order flow, we also assume that the libration frequency is larger than twice 
the mean rotation rate (spin frequency). This hypothesis guarantees that no inertial 
waves are excited and no critical latitude singularity is present. It will allow us to obtain 
an expression for the mean zonal flow in the whole shell. An expression for the mean 
zonal flow will also be provided for smaller frequencies (smaller than twice the mean 
rotation frequency) but it will exhibit a divergence associated with the presence of a 
critical latitude. We shall see that for any libration frequency, the presence of an inner 
core significantly modifies the mean zonal flow and generates so-called Stewartson layers 
( Stewartson 1966) across the cylinder tangent to the inner core. 

The paper is organized as follows. In section 2, we present the mathematical framework 
of the problem. We provide the governing equations and the hypothesis for the asymptotic 
analysis. The leading order solution at the libration amplitude e is derived in section 3. 
In section 4, we focus on the mean flow correction at the order e 2 . We explain how the 
nonlinear interactions within the Ekman layer generate a mean zonal flow in the bulk 
of the shell. We show that a discontinuity appears across the cylinder tangent to the 
inner core which re quires the introdu ction of viscous layers as for differentially rotating 
concentric spheres ( Stewartson! [l966). The structure of the solution in these layers is 
provided. We show that both axial vorticity and axial flow are the strongest in these 
layers. The details of the calculation leading to the expression of the different fields are 
given in the appendix. In section 5, the asymptotic solution is compared to numerical 
results obtained for both a sphere and a shell and a fairly good agreement is obtained. 
A brief conclusion with some perspectives is provided in section 6. 



2. Framework of the asymptotic analysis 

We consider a homogeneous and incompressible fluid of kinematic viscosity v contained 
in a spherical shell of external radius R ex t and inner radius Ri n t ■ This shell rotates in the 
laboratory frame at angular velocity which sinusoidally oscillates around a mean angular 
value flo- We use R ext and flo^ 1 as the length scale and the time scale respectively. As 
both the libration forcing and the geometry are axisymmetric, we shall assume that the 
flow remains axisymmetric, that is its velocity and pressure field are independent of the 
azimuthal variable (p. It is then convenient to use the spherical coordinates (r, 9, </>) and 
to express the velocity field (it r , ug, u^) as functions of tp and x defined by: 



f 



dtp 

" 2 sin 9 89 ' 



f 



Ug 



dtp 

r sin 9 dr ' 



X 



r sin 9 



(2.1) 



The use of a toroidal/poloid al decomposition of the velocity field is al so possible (see for 
instance Dormv et al. 1998 ). but we have preferred the formulation of Stewartson] ( 1966h 
to use the similarity with his analysis. 

In the frame rotating at the mean angular velocity, the dimensionless Navier-Stokes 
and continuity equations can be expressed as 



dD 2 tP 
dt 



cos 9 — — E D tp = 
dr 



r d9 
dtp dD 2 iP 



2 x ( sin 9 dx 



r 2 sva.9 \ dr dO 



dtp dD 2 tP 
~d9 dr 



r 2 sin 2 i 



2D 2 %P 
r 2 sin 2 6 



r d9 
sin e dtp 



dr 



r 



de 



dtp 
cose—- 

dr 



(2.2) 
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(2.3) 
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where 



o d 2 sin 9 d ( 1 d , 
o at* V sin cw ; 



and -E is the Ekman number defined by: 



v 



Rext 2 

We consider no-slip boundary conditions for a general configuration where the outer 
sphere and the inner core do not librate at the same amplitude. The boundary conditions 
in the rotating frame of reference then read for the functions tp and x- 

dip 

— = ip = and X = e sin 2 cos(cjt) at r = 1 (2.6a) 
— — = tp = and x = £ a fl2 sin 2 $ cos(wi) at r = a (2.66) 

where a = Ri n t/ Rext is the aspect ratio of the shell, e the amplitude of libration of the 
outer sphere and uj the frequency of libration. The parameter a measures the relative 
amplitude of libration of the inner core. For a non-librating inner core, we impose a = 0, 
whereas when it is librating at the same amplitude as the outer sphere, a = 1. For the 
case of a sphere, without inner core, only the condition at r = 1 has to be applied. 

The objective of our work is to characterize the steady flow induced by the librational 
forcing. We assume that the forcing amplitude s is small such that the flow remains 
stable with respect to the centrifugal instability. Moreover, we consider the limit e<l 
in order to use asymptotic methods. In this limit, we can expand the functions ip and \ 
in power of e: 

ip = sip 1 +e 2 ip 2 + o(e 2 ), (2.7a) 
X = e X i+e 2 X2 + o(s 2 ). (2.76) 

The first order solution is expected to be the linear response at the librating frequency 
forced by the boundary while the second order represents the solution generated by the 
nonlinear interactions of the first order solution with itself. As in all weakly nonlinear 
analysis, both harmonic corrections oscillating at 2w and mean flow corrections (the zonal 
flow) are expected. It is this second part owing to its peculiar structure that will be our 
point of focus. 



3. Linear solution 

In this section, we consider the solution obtained at the order s. We are interested in the 
solution for very small values of E, as encountered in geophysical applications. For this 
reason, it is natural to consider E as an asymptotically small parameter. Nevertheless, 
this limit is taken after having considered e — > 0. In the small E limit, we expect viscous 
effects to be localized in space. In our analysis, we shall require that, for the first order 
solution, they are localized near the boundaries. This strong constraint can be satisfied 
if no inertial waves are excited in the bulk by the harmonic forcing. This is verified if the 
frequency of the forcing is outside the range of inertial wave frequencies, that is uj > 2, 
which is the condition that we shall assume in the following. 

For uj > 2, the linear response is expected to be mainly localized close to the bound- 
aries. This solution can be obtained by performing a classical asymptotic analysis as 
done for Stokes or Ekman layers. Let us consider first the boundary layer on the outer 
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sphere at r = 1 . The boundary layer solution is obtained by introducing the local variable 
r = (1 — r)l\f~E. Let us write the leading order expression for if) and \ as 



Xi 



E^[ o) e wt + c.c. 

(o) iujt , 

Xi e + c.c. 



(3.1a) 
(3.16) 



where the superscript (o) indicates that we consider the function in the viscous layer at 
the boundary of the outer sphere and c.c. denotes complex conjugate. Inserting (I2.7a .b) 
with (I3.la .b~) in (12.21) and (12.3[) . we obtain, at leading order in E and e, the equations 



dr 2 

(o) 



2 cos9 
- 2 cos 9 





d A 4 o) 


dr 


dr 4 




d 2 x ( ° ] 


dr 


dr 2 



with the no-slip boundary conditions 

4 o) (r = 0,8) = , 
d f 4°\r = 0,0) = , 
x ( o) (f = 0,6») = (sin 2 6>)/2 . 

These equations together with the boundary conditions give 

i sin 2 9 II -cr x + f 1 - e" A - 



;(°) 
^1 



(o) 

x\ 



X 
sin' 



-A+ r 



-A- f 1 



(3.2 a) 
(3.26) 



(3.3a) 
(3.36) 
(3.3c) 



(3.4 a) 
(3.46) 



where we have defined 



(1+i) 



± cos 9. 



(3.5) 



This expression can also be found in Bussei (1968). Note that the condition uj > 2 
guarantees that the boundary layer solution does not diverge for any value of 9. Such 
a divergence, called critical latitude singularity, is known to be a source of inner shear 
layer s as documented for the case of precession (see for instance iHollerbach fc Kerswell 
19951) . 



A similar solution is obtained in the viscous layer near the inner sphere using the local 
variable r = (r — a)/ \f~E~: 



/(») 



9 . 9 

ia sin 



Xi 



-A+r 



-A_ r 



a 2 sin 2 



re- A +' 



(3.6a) 
(3.66) 



As f (resp. f) goes to infinity, Xi (resp. x*' '') g° es to zero while ip^ (resp. -0^) goes 
to a non-zero constant. This part of the solution is responsible for a contribution in the 
bulk at the order e VE which corresponds to the Ekman pumping. 
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4. Mean zonal flow 

4.1. Flow asymptotic structure 

We now consider the nonlinear solution at order e 2 . This second order correction is 
created by the nonlinear interaction of the first order solution with itself. Note first that 
the nonlinear interaction of the solution in the bulk (the Ekman pumping) is expected 
to generate a correction at order e 2 E. As we shall see, this order is actually smaller that 
the order of the mean flow correction that we will find in the limit E — > 0. The main 
contribution in the bulk will be found to come from the nonlinear interactions in the 
boundary layers and to be of the form: 



iP2~VEi> 2 (r,6) (4.1a) 

X2~X 2 (r,0) (4.16) 

where ip2 and X2 are both independent of time and E. 

Equations (j2.2[) and (|2.3[) imply that ip 2 and X2 would satisfy in the bulk 

d sa\9 d \ - „ , / „ 8 sin# d \ „ . , „. 

>2=0 and [cose- — U 2 = 0, (4.2) 



Or r 89 J \ dr r 89 

which can be rewritten in cylindrico-polar coordinates (p, <f>, z): 

^1 =0 and ^=0. (4.3) 

8z 8z 

This means that, in the bulk, these two functions would only be dependent of the radial 
coordinate p — r sin 9 and could be written as 

&(r,0)=¥ a (p), (4.4a) 
X2(r,e)=X 2 (p). (4.46) 

This property is nothing but a consequence of the Taylor-Proudman theorem. 

In the next section, we calculate the functions ^(p) and X%(p) by performing the 
matched asymptotic analysis of the solutions in the boundary layers and in the bulk. The 
asymptotic structure of the problem is provided in figure [1] In addition to the boundary 
layers, the volume is split in two main regions I and II. Because one of the boundary layer 
disappears when we cross the cylinder tangent to the inner core, different approximations 
are obtained in regions I and II. We shall see that the functions are discontinuous across 
the boundary between regions I and II and that smoothing this singularity requires thin 
nested layers (region III) in which viscous effects have to be introduced. 

In the following, we only consider mean flow quantities: the bars are then dropped out 
(e.g. xi w iH be noted xa)- 

4.2. Boundary layer solution near the outer boundary r = 1 

Using the boundary layer variable f = (1 — r)/VE, the equations (|2.2|) and (|2.3[) for the 
stationary component in the boundary layer at order e 2 write 

^2_^2cos0^|-=iv4 o) , (4.5) 



8f A df 
or z or 



_ 2 cos 9 = M4 o) , (4.6) 




Figure 1. Structure of the flow in the librating shell showing the presence of three different 

zones. 



where NL[ o) and NL%> are the stationary part of the nonlinear forcing terms given by: 



-(c) 



NL^=- 2C ° S9 



sin 



d 2 tp[ 0> dip] ' (o) d\ 



df 2 df ' 1 df 
d4 o) d X [ o) d^ { ° ] d X [ o) 



NL 



(o) 



sim 



df dOdf 2 86 df 3 



df d9 



89 df 



From (|4.5|) and (|4.6|) . we obtain a fifth-order equation for i/^ 



(°). 



8_ 

df 



d*_ 

df 4 



4 cos 



( o) = dNL^_ _ 2 cos0M7 (o) = ^(o) 
<9f 



This equation has to be solved with the boundary conditions: 

4 o) (r = 0) = 0, M-(f = 0)=0, lim Va o) (*0=*a(Bin( 

C/r r->+oo 
» 



(4.7) 
(4.8) 

(4.9) 
(4.10) 



The expressions of M> } and NL\°' are cumbersome. We can write them in a compact 
form as 



NL^=J2 [A^ + fB^)] 
i=i 

NLP^idW+fDtiO)] 



-Ml T 



-fil r 



C.C. 



1 = 1 



where the exponents /i/ are related to A± defined in (|3.5[) by 
jUi = A+ + A*, /i 2 = A + + A"i_, fi3=X- + \*_, H4 = X- 



M5 = A_ 



(4.11a) 
(4.116) 

(4.12) 



Expressions for the coefficients Ai(9), Bi{9) , Ci(9) and Di(9) are provided in appendix A. 



A. Sauret & S. Le Dizes 
The solution to (|4.9|) which satisfies the boundary conditions (|4.10p reads 



V4 o) (f,6>) = 4 (M) + *2(sin< 



1 



(4.13) 



where 



E 



4 cos 2 9 + m A 



Bi + mAi 



4 COS 2 + 



-Mi'" 



5/ 



^(4 COS 2 6» + /i; 4 ) 



+ c.c. 



and where we have defined 



From (14.51) , we obtain x^ > which must satisfy the matching condition: 

lim x { 2° ) (r)=X 2 (sm9) 1 

r— >-\-oo 



k = (1 + i)vcosl 



(4.14) 
(4.15) 

(4.16) 



as 



x 2 o) (^) = 



1 



2 cosfl 



5 3 V. (o) 




+ c.c. 



The above expression satisfies the no-slip boundary condition x^ (0, 0) 
and A" 2 satisfies the relation: 



X 2 (sm9) = T(sin9;ui) - 2Vcos 9 * 2 (sin t 



where 



J-(siii9]Uj) 



1 



2 cos( 



d 3 (j) 



df 3 




Pi + m Q 
m 2 



+ c.c. 



-A? 2 (sin0). 

(4.17) 
= only if \I> 2 

(4.18) 



(4.19) 



Equation (|4.18l) is a constraint obtained from the matching between the solution in 
the bulk and the solution in the boundary layer on the outer sphere. 

A similar constraint is expected from the boundary layer on the inner sphere. The 
same analysis in this boundary layer leads to 



(4.20) 



X 2 (asin9) = a 2 a 2 „F(sin 9: ui) + 2Vcos1? ^(asinfl). 

4.3. Solution in the bulk 

Relations (|4.18p and (|4.20l) permit us to obtain the functions ^ and X 2 which charac- 
terize the mean flow correction at leading order in the bulk. 

We first consider region I, defined by p > a (see figure [T]). In this region, the only 
contribution comes from the viscous layer close to the outer sphere. However the flow in 
this region must satisfy a condition of symmetry on the equator: the axial flow should be 
antisymmetric with respect to the equator and therefore should vanish on the equator. 
Since it does not depend on the axial coordinate, this condition implies that it should 
vanish everywhere in region I. In other words 

* 2 O0=O, (4.21) 

from which we deduce from (|4.18j) . that for p > a (i.e. in region I): 

X 2 (p)=F(p;u). (4.22) 
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In region II, both the conditions of matching with the inner and outer boundary layers 
have to be considered. Replacing s\\\9 by p in (|4. 18|) and asin# by p in (|4.20|1 . lead to 

X 2 (p) = -2 [1 - p 2 ] 1/4 tf 2 (p) + F{p- u), (4.23a) 

X 2 (p) - 2 [1 - (p/a) 2 ] 1/4 * 2 (p) + a 2 a 2 ^(p/a; w), (4.236) 

from which we deduce, for p < a (in region II) 

= ^(p;^)- a W(p/ a ;^) 

2(l-p 2 )V4 + 2[l-(p/a)2] 1 /4' ^ ' 

, N a 2 a 2 (l-p 2 ) 1 /4j-(p/ a;w )+ [l-(p/a) 2 l 1/4 J-(p; W ) 

(l-p 2 )l/4 + [l_( p / a )2]l/4 

Using (12.11) . we can also calculate from 9 2 and X 2 , the mean azimuthal velocity u^, the 
mean axial velocity u Z2 and the mean axial vorticity uj Z2 which satisfy in regions I and II 

X 2 VEd^2 ldX 2 

U <t>2 = 1 U *2 = ) W 22 = ( 4 ' 25 ) 

p p op p op 



Note that expressions (|4.24b .b) resemble the expressions obtained bv lProudmanl ( 1956h 



for the flow between two spheres rotating at different speeds. The difference is in the nu- 
merator which is now a more complicated function associated with the nonlin ear interac- 



tion in the boundary layer of the flow generated by libration. However, as in iProudman 



1 1956t ). V&2 and X 2 exhibit a discontinuity across p = a. In particular, note that 

X 2 {p) ~ X+=T{a-uj), (4.26a) 



X 2 (p) ~ X- = a 2 a 2 .F(l;w) = — — . (4.266) 

2 oj z 

The singulariti es of ^2 and X 2 can be smoothed across a series of viscous layers as 
first shown by IStewartso n (1966). In the next section, we provide approximations of the 



solution in these layers. 

It is worth noting that the solution in region II can be deduced for the solution in region 
I. In other words, if we exclude the Stewartson layers, the solution in a spherical shell 
can be deduced from the solution in a sphere. Indeed, all the information is contained 
in the function T{j)\bS) which defines X 2 in region I (or in a sphere). Unfortunately, we 
have not been able to obtain a simple expression for this function [see formula (|B ip given 
in appendix B]. In figure [2 we have plotted T(p;ui) versus p for various values of the 
libration frequency oj. A contour plot of J"(p; w) in the (p, ui) plane is also shown in figure 
[3] Interestingly, we observe that T changes sign as p varies when u> < 8.5. This means 
that the mean azimuthal flow generated by libration could be anticyclonic close to the 
axis and cyclonic further away. 

For large libration frequency oj, an asymptotic expression for J-(p;uj) is obtained as 



T( \ P 2 V^ y/2p 2 (l-p 2 ) 3 / 4 , p 2 (7p 2 -5) /1\ 

j- (p;w) _ __ + 4w2 +o[^ (4.27) 

uj— too 

from which we can easily deduce the behavior of X 2 , u<j, 2 , ^2 and u Z2 . In figure EJ we can 
observe that (|4.27p provides a very good approximation for T for ui = 10 and larger. 
Close to p = 0, the function J"(p;w) can be related to the expression obtained for a 
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0.14 




0.2 0.4 0.6 0.8 1 

P 



Figure 2. The function J-(p; uj) versus p for different values of the libration frequency u). This 
function characterizes ^2 and X% in the bulk. Expression (|4.27[) is plotted in dash gray for 
uj = 10 and uj — 20. 




P 

Figure 3. Contour levels of the function !F{p\uj) in the (p, uj) plane. 



1970]). We have T(p;ui) ~ tt C yi{u})p 2 as p — > which leads to 

Xi ~ { -^^n Cyl (u:)p 2 . (4.28) 

p->0 

The function flc y i (w) corresponds to the mean angu lar velo city generated by the libration 
of a cylinder. Its expression was first obtained by IWane 
2012t) . Note the simple expression for the function T at p = 1: T(1;lu) = 1/(4cj 2 ). 
It is also interesting to provide expressions for ^2 and X% in the limit of thin shell 



1970) (see also ISauret et al 
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Figure 4. Structure of the Stewartson layers. The inner layer has a thickness of E 1 ^ 3 and the 
external and internal outer layers have a thickness of E 1 ^ and E 2 ^ 7 respectively. 



(a 



*9 



(1 



'4(1-^2)1/4 



+ 0(1 -a), 



a 



(4.29a) 



(4.296) 



Expression (|4.29b ) implies that ^2 becomes small of order (1-a) when inner and outer 
sphere are librating at the same amplitude (a 2 = 1). Expression (|4.29b ) means that the 
p dependence of the function X2 in the thin shell limit is the same as for a full sphere. 

4.4. Solution in the Stewartson layers 



The structure of the viscous layers is the same as in IStewartsonl (|1966i ) and is sketched 
in figure 2J There arc three different layers around p = a, two outer layers scaling as 
\p — a\ = 0(E 1 / 4 ) for p > a and \p — oj = Q(E 2 / 7 ) for p < a, and an inner layer 



where \p — a\ = 0{E 1 / 3 '). As shown by IStewartson ( 1966 ). the outer layers guarantee 
the continuity of \i (and u,p 2 ) and of its first derivative, while the inner layer is for the 
smoothing of ip2 (and u Z2 ) and higher derivatives of X2- More details on the calculation 
leading to approximations of the solution in each layer are given in appendix C. In this 
section, we only provide the leading approximation in each layer. In the appendix, we 
explain how the next order which is 0{E 1 / 2S ') can be obtained. 

In the outer layers, simple expressions can be obtained. In the external outer layer 
(p > a), \2 varies exponentially between the two values X ± = lim p _;. a ± X^{p) of X2 on 
either side of a in regions I and II: 



(l- e -P)+X-e-P, 



with 



P=(p~a)/[(l-ay*E^ 
In the internal outer layer (p < a) , it is a constant 



(4.30) 
(4.31) 

(4.32) 



The function X2 (and therefore 11^2) is thus of same order in the Stewartson layers than 
in the bulk. By contrast, the vorticity associated with the mean flow correction is larger 
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Figure 5. Axial vorticity uj Z2 /u} Z2 (a) as a function of p = (a - p) [a/(2 (1 - a 2 ) 2 E 2 )] 1/7 for 
p < a (internal outer layer), and p = (p — a)/ (1 — a 2 ) 3 / 8 _E 1//4 for p > a (external outer layer). 



in the outer layers. Its main component is the axial component uj Z2 which reads 

u z ~ _E _1 / 4 — , e~ p in the external outer layer, (4.33a) 

a(l — ar) 1 */ 8 

oj Z2 ~ iT 1/4 — —-7- — — in the internal outer layer, (4.336) 

a (1 — a^' 8 1 (0) 

where the internal outer layer variable p is defined by 

"I 1/7 



(4.34) 



The function f is obtained in the appendix. It can be expressed in terms of Bessel 
functions [see expression (|C 16|) ]. In figure [SJ we have plotted w Z2 /w Z2 (a) in each outer 
layer. Note that the vorticity is strongly peaked at p = a where it reaches its maximum 
value (in amplitude) 

"-<"> = E ~ v, i^?w ( " 5) 

The discontinuity of the derivative at this point is smoothed in the inner layer. 

Contrarily to the azimuthal velocity and the axial vorticity, the function ?/>2 and the 
axial velocity depend on the axial variable z at leading order. This dependence is linear 
in the outer layers [see expressions (|C 8[) and (jC 18|) in the appendix C] but it is more 
complex in the inner layer. It is in this inner layer where \p — a\ = 0{E X ^) that ^2 
and u Z2 are the largest. They scale respectively as -02 = 0(E 19 / 42 ) and u Z2 = 0(E 5 / 42 ). 
Note in particular that the axial flow is much larger than in region II where it is 0(^/E). 
Unfortunately, the expressions for ip2 and u Z2 are not as simple as for \2- In the inner 
layer (\p — a\ = 0(E 1 ^ 3 )) : ip2 and u Z2 can be written as 

z 
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(a) (b) 

Figure 6. The functions (a) \{Z, R) and (b) J(Z, R) versus R for different values of Z 
(Z — 0.1 : 0.1 : 0.9). These two functions provide the spatial variation in the inner layer, 
of tp2 and u Z2 respectively. 




where I and J are the functions given in the appendix by (jC 23|) and (jC 28|) and A 3 and 
A 4 are related to SX = X + - X~ = F(a;uj) - a 2 a 2 /(Auj 2 ) by 

a a 3 / 28 

A 4 (a,o;) = a25/28 (1 ^ fl2)25/84 5X(a,u), (4.37b) 

with a a = -l/[2 6 / 7 f'(0)] w 0.45. 

The variations of the functions I(Z, R) and 3(Z, R) with respect to R at various height 
Z are shown in figure|ni Note that both functions are the largest for small values of Z, that 
is close to the inner sphere. It is also interesting to note the peculiar form of the axial flow 
profile, which is composed of a main jet surrounded by two weaker counter-propagating 
jets. 

It is important to point out that tp2, u Z2 and w 22 depend on w and a only via the 
function 5X and that they are all proportional to that function. When the inner core is 
not librating, a = 0, and 5X reduces to the function T(a;ui) which has been plotted in 
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figure [3l The libration of the inner sphere reduces the contribution from the outer sphere 
by a simple quantity a 2 a 2 /(4w 2 ). If only the inner core was librating (the outer sphere 
not librating) and if we had used for e the amplitude of libration of the inner sphere, we 
would have only obtained the inner sphere contribution 8X = — a 2 /(4w 2 ). The contours 
of this function in the (a,w) plane are plotted in figure [71a). This plot which describes 
the contribution from the inner sphere has to be compared with the one obtained for the 
outer sphere contribution shown in figure [3] We clearly see that the contribution from 
the inner sphere is always of the same sign whatever a and oj: therefore, the libration 
of the inner core always generates a negative (anticyclonic) mean axial vorticity and a 
main axial flow which goes toward the inner sphere. In figure [T^b), we have plotted the 
variations of 6 X with respect to both variables a and u) when inner and outer spheres 
librate at the same amplitude (a = 1). Note that if we compare this plot with figure El 
we observe that the inner sphere provides a non-negligible contribution as soon as the 
radius a of the inner core is not small. In particular, it strongly increases the domain of 
parameters where the mean axial vorticity is anticyclonic. 



5. Numerical comparison 

In this section, the results of the asymptotic theory are compared to nu merical results 
obtain ed for finite values of e and E. Both published results obtained bv lCalkins et all 
(2010) and results obtained using the finite elements commercial software Comsol Multiphysics® 
are considered. Both codes rely on a 2D (axisymmetric) model to be able to reach E as 
small as E = 10~ 6 . For more details about our own numerical procedure, the reader is 
referred to Sauret et all (|2010l |20 12h where the numerical procedure has already been 
used for similar problems. 

The first and easiest comparison that can be made is for the libration in a sphere. 
In that case, there are no Stewartson layers and no axial flow in the bulk. The main 
contribution to the mean flow correction is an azimuthal flow which is theoretically 
expected to be 

u^e*^l, (5.1) 
P 

for small Ekman number E and small libration amplitude e. 

Thus, the theory predicts that the zonal flow should scale as the square of the libration 
amplitude and be independent of the Ekman number. These properties have been tested 
by comparing numerical results obtained with Comsol Multiphysics to the theory for e 
ranging from 0.01 to 0.3 and E ranging from 10~ 3 to 4.10~ 5 . 

In figure [HIa), we display the results obtained for fixed frequency and amplitude and 
different Ekman numbers. We can observe that the numerical curve which is already close 
to the theoretical curve for E = 10 -3 gets even closer as the Ekman number decreases. 
The departure obtained close to p = 1 is associated with the boundary layer on the outer 
sphere which is not included in the theoretical expression (|5 . 1 1) for the solution in the 
bulk. When E decreases, the boundary layer which scales as E 1 / 2 becomes thiner, as 
observed. The amplitude e has also been varied but its effect has been found very weak. 
For instance, the curves obtained for e = 0.01 and e = 0.3 cannot be distinguished from 
those plotted for e = 0.1 in figure a). 

In figure [UJb), we have compared the theory to numerical results for different libration 
frequencies. Again, we observe a very good agreement for all frequencies. These results 
provide a validation of the asymptotic analysis that has led to expression (|5.1[) . 
In a librating shell, we expect from the theory that the azimuthal zonal flow remains 
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' 0.2 0.4 0.6 0.8 1 

Figure 8. Mean azimuthal flow u ( j >2 — u^e~ 2 in a sphere. Comparison of numerical results 
(thin lines) obtained in the equatorial plane (z = 0) with theoretical results (thick line), (a) For 
fixed frequency and amplitude (here cj = 3; e = 0.1) and different Ekman numbers (Dot-dashed 
line: E = 10" 3 ; Dotted line: E = 5.10" 4 ; (Thin) solid line: E = 10" 4 ; Dashed line: E = 4.10" 5 ); 
(b) for fixed Ekman number and amplitude (E = 5.10 -5 , e = 0.01) and different frequencies. 



independent of z in the bulk but exhibits a discontinuous behavior across the cylinder 
tangent to the inner sphere at p = a. Such a discontinui ty is visible in the n umerical 
results. In figure HJa) , the numerical results obtained by ICalkins et al. ( 201fj| ) for two 
values of the Ekman number (E — 10~ 5 and 1CP 6 ) are compared with the theoretical 
prediction in the bulk (regions I and II). 

The agreement is good if we do not consider the region close to p = a = 0.35 where the 
theoretical solution in the bulk is discontinuous. In particular, note that for both Ekman 
numbers, and any value of z, the azimuthal velocity profiles collide far from p = 0.35 on 
a single curve, as expected from the theory. As in figure [51 the fast localized variations 
observed on the numerical curves correspond to the boundary layer on the outer sphere 
p 2 + z 2 = 1. Their position changes because different z are considered. In the Stewartson 
layers, the solution becomes dependent of the Ekman number. In figure[9jb) the numerical 
results for z — 0.6 are plotted with respect to the outer layer variables (p — a)/E 1 / 4 and 
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Figure 9. Mean azimuthal velocity for e = 2.2 10~ 3 , a = 0.35, u = 2.2 and a = (no libration 
ot the inner core), (a) Solution in the bulk: u ( p 2 = u</,£ -2 versus p. The thick solid line is the 
theoretical prediction in the bu lk (regions I and II). The thin lines are the numerical results 
obtained bv lCalkins et al\ (|201Gi ) at three different z (z = 0.4, 0.6, 0.8) and two different Ekman 
numbers (E — 10 -5 : dashed lines, E — 10 -6 : solid lines), (b) Solution in the Stewartson layers: 
u r f >2 = u</,e -2 versus (p — a)/_E 1//4 for p > a and versus (p — a)/E 2 ^ 7 for p < a. The thick solid 
line is the leading order theoretical prediction in the outer Stewartson layers. The thick dashed 
line is the theoretical prediction for E = 10 -6 where the correction in E 1 ^ 28 has been taken into 
account. The thin lines are the numerical results of ICalkins et al\ (|2010T ) at z — 0.6; solid line: 
E = 10" 6 , dashed line: E = 10" 5 . 



(p — a)/E 2 / 7 . These results are compared to the leading order theoretical prediction as 
well as to the two-order approximations obtained for E = 10~ 6 in this plot. We can first 
note that there is a noticeable effect of the second order correction on the theoretical 
curve. This correction which is of order E 1 / 28 is clearly not small for an Ekman number 
E = 10~ 6 . If this correction is taken into account, a relatively good agreement is obtained 
between the theoretical curve and the numerical results close to p = a. 
A similar comparison is performed for the axial vorticity in figure [10] and for the axial 
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Figure 10. Mean axial vorticity for e = 2.2 10~ 3 , a = 0.35, u> = 2.2 and a = (no libration of 
the inner core), (a) Solution in the bulk: uj z e~ 2 versus p. The thick solid line is the theoretical 
prediction in t he bu lk (regions I and II). The thin lines are the numerical results obtained by 
Cal kins el al\ i|2010f ) at three different z (z = 0.4, 0.6, 0.8) and two different Ekman numbers 
(E — 10~ 5 : dashed lines, E = 10 -6 : solid lines), (b) Solution in the Stewartson layers: uj z e~ 2 E 1 '' 4 
versus (p — a)/E x ^ A for p > a and versus (p — a)/E 2 ^ 7 for p < a. The thick solid line is the 
leading order theoretical prediction in the outer Stewartson layers. The thick dashed line is the 
theoretical prediction for E = 10 -6 where the correction in j? 1 / 28 has been taken into account. 
The thin lines are the numerical results of lCalkins et al\ (|2010T ) at z = 0.6; solid line: E = 10~ 6 , 
dashed line: E = 10" 5 . 



velocity in figure [TT] In figures [TDT a) and lllf a). we have compared the numerical results 
with the theoretical solution in the bulk, that is (14.251) with (|4.24b .b) in region II (r < a), 
and with (|4.21|) . (|4.22l) in region I (r > a). A good agreement is observed far from p = a. 
Note in particular that both axial vorticity and axial velocity are independent of z far 
from a as expected from the theory. However, there are systematic differences which 
were not present in the sphere. We suspect that these differences are due to 0(E 1 ^ 28 ) 
corrections generated in the Stewartson layers, and which are not negligible for E = 10~ 6 . 
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Figure 11. Mean axial velocity for e = 2.2 10 , a = 0.35, w = 2.2 and a = (no libration 
of the inner core), (a) Solution in the bulk: u z e~ 2 E~ 1 ^ 2 versus p. The thick solid line is the 
theoretical prediction in the bu lk (regions I and II). The thin lines are the numerical results 
obtained by Calkins et al.\ (|2010T ) at three different z (z = 0.4, 0.6, 0.8) and two different Ekman 
numbers (E — 10 -5 : dashed lines, E — 10 -6 : solid lines), (b) Solution in the inner Stewartson 



layer at 



0.7: 



"E 



-5/42 



us (p — a)/£ ,1//3 . The thick solid line is the leading order 
theoretical prediction in the inner Stewartson layer. The thick dashed line is the theoretical 
prediction for E = 10 -6 where t he correction in jg 1 / 28 has been taken into account. The thin 



lines are the numerical results of ICalkins et al\ (|2010T ) at z = 0.7; solid line: E = 10 
line: E = 10" 5 . 



dashed 



In figure fTDTb). the axial vorticity is plotted with respect to the outer-layer variables. 
Although a slight shift is observed in the numerical results, we can see that the vorticity 
peak with its scaling in _E -1 / 4 is qualitatively described. In ngurcfTTTb). the axial velocity 
is plotted with respect to the inner layer variable. The agreement of the numerical results 
with the second-order theoretical approximation is remarkable for this value of z. Again 
a slight shift is observed which seems to diminish as E decreases. Other values of z have 
also been tested and a less good agreement has been observed for smaller values of z. 



Libration-induced mean flow in a spherical shell 



19 



The numerical results shows a weaker jet whereas the theory predicts the opposite. We 
suspect that this qualitative difference could be due not only to the large value of the 
Ekman number but also to the particular structure of the Stewartson layer for small z 
which has not been resolved here. 



6. Conclusion 

In this paper, we have analyzed the mean flow induced by longitudinal libration (that 
is small oscillation of the rotation rate) in a spherical shell and in a sphere, when the 
oscillation frequency is larger than twice the mean rotation rate. Using asymptotic meth- 
ods in the limit of small Ekman numbers E and small oscillation amplitude e, we have 
been able to show that the main contribution to the mean flow is, in the bulk, an az- 
imuthal flow which scales as e 2 , remains independent of the Ekman number and depends 
on the cylindrical variable p only. We have shown how such a zonal flow is generated by 
nonlinear interaction of the harmonic solution in the boundary layers. The dependence 
of the zonal flow with respect to the libration frequency oj and the ratio of the core size 
a has also been comprehensively analysed. In particular, we have shown how the zonal 
flow in a shell can be deduced from the zonal flow in a sphere. We have also seen that 
for moderate frequencies 2 < u < 8.5, both cyclonic and anticyclonic mean rotation can 
be created by libration, whereas only weak cyclonic rotation is created for uj > 8.5. 

The presence of the inner sphere has also been shown to create a discontinuity of the 
azimuthal flow across the tangent cylinder at p = a which is smoothed in a series of nested 
viscous layers of widths E 1 / 4 , E 2 / 7 and E 1 / 3 , as for the configuration of differentially 
rotating spheres ( Stewar tson 1966). The mean flow correction has been calculated in each 
of these layers. We have obtained that an axial vorticity field of order e 2 E~ x I 4 is created 
in the outer Stewartson layers, and that a non-uniform axial flow of order e 2 E 5 / 42 is 
present in the inner layer (|p — o| = 0(E 1 ^ 3 )). We have shown that the spatial structures 
of these two fields in the Stewartson layers are universal functions which do not depend 
(after rescaling) on any parameter. The ratio of the core size, the oscillation frequency 
as well as the relative oscillation amplitude of the inner sphere, have been shown only 
to intervene in an amplitude factor which has been fully characterized. In particular, we 
have demonstrated that the libration of the inner sphere adds a very simple quantity 
to the amplitude factor obtained without inner sphere libration, which always tends to 
increase the anticyclonicity of the vorticity field, and to add an axial flow from the inner 
core to the outer sphere. 

The asymptotic r esults have been compared to numerical results of the literature 
(jCalkins et a/.ll2010l ) or obtained with a finite-element code, and a good agreement has 
been demonstrated for the zonal flow in the bulk. A reasonable agreement has also been 
obtained for the peculiar structure of the solution in the Stewartson layers despite the 
largeness of the Ekman number in the simulations. In particular, we have shown that 
it was necessary to consider the second-order 0(E 1 ^ 28 ) correction to obtain a correct 
agreement. It is worth mentioning that our theory is valid up to 0(E 1 ^ 14 ). Even for 
realistic astrophysical values of E {E » 10~ 14 ), we cannot guarantee that the error will 
remain small. This constitutes a weak point of the theory that is important to keep in 
mind when the results are applied. 



7. Discussion 

The present paper has considered oscillation frequencies larger than twice the rotation 
rate in a spherical shell. We can naturally address the question of the libration response 
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Figure 12. The function J-(p;uj)/p 2 versus p for different values of uj between and 2 . The 
function T diverges at p c = yl — uj 2 jl- Expressions (|B II) and (|B3|) have been us ed for T for 
uj > 0. The dashed line for u; = is the expression (|7.2I) obtained by IBussel (|2010al ). 



when the frequency is smaller than twice the rotation rate or when the boundaries of the 

shell are elliptically deformed (that is the outer/inner sp heres are ell i psoids). 

The first point has already been addressed in a c ylinder (IWandll97filSauret et aZ.ll2012h 
and for vanishingly small frequencies in a sphere (|BusseH2010ct ISauret et al. I I2010D . For 



this point, the first difficulty co mes from the fact that inertial m odes can be resonantly 
excited by the harmonic forcing (jMorize et al. 2010l ; TilgnerlfeOQ? ). When this occurs, the 
zonal flow is strongly modified and is no longer dominated by the nonlinear interaction in 
the boundary layers (jSauret et aLll2012n . The second difficulty is the singular behavior of 
the boundary layer solution when twice the cosine of the inclination angle of the boundary 
matches the oscillatio n frequency. This singular ity is known to generate internal shear 
layers in the bulk (see lKerswelllll995t lKidall2011l) whose effects on the mean flow remain 
to be quantified. 

If we neglect the effect of such internal shear layers and assume no resonance, the mean 
zonal flow can be calculated using the same approach as for u > 2. The expressions (|4.21l 
I4.22p and (|4.24b .b) for <F 2 and X 2 still apply but expression (|B Q for F(p; uj) obtained 

uj 2 /A. For p < p c 



for uj > 2 can only be used for p > p c = yl — uj 2 /4. For p < p c , expression (IB 31) has to 
be used. This expression is obtained by the same analysis by changing the definition of 
the square root in ()3.5j) . Both expressions diverges at p c which corresponds to the radius 
of the critical latitude on the outer sphere. In a shell, this singularity in J- implies two 
singularities in X2. The second one which is associated with the critical latitude on the 
inner sphere is at p c = ay/1 — uj 2 /4. From an asymptotic point of view, these singularities 
mean that another scaling has to be found for X2 close to each singularity. 

In figure Q21 we have plotted T/p 2 , that is the angular velocity of the zonal flow 
in region I (see figure [TJ, for different values of uj. It is interesting to note that for 
frequencies ranging from to 1, we obtain almost the same constant retrograde rotation 
for all p < 0.6. 

As uj —> 0, expression (|B 3j) reduces to 



(59p 2 -72)p 2 
HP ^ ] ~ 480 (1 -p 2 ) • 



(7.1) 



This 



expression slightly differs from the expression obtained bv IBussel ( 2010a ) which 
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reads 

(51.8 P 2 -72)p 2 
HP ^ ] ~ 480 (1 -p 2 ) • (7 ' 2) 

Expression (|7.ip is valid when e <C w <C 1 whereas Busse expression applies when lj <C e. 
As noticed in figure I12[ the difference between these two expressions is particularly 
visible for p > 0.7. 

The small frequency limit is interesting because small frequency inertial modes are 
generally too damped to be resonantly excited. Moreover, for small lj, internal shear 
layers are limited to the very small region near the equator and to the neighborhood of 
the cylinder tangent to the inner sphere, that is within the Stewartson layers. Therefore, 
they are disconnected from the bulk regions I and II where our solution is expected to 
apply. 

The second po int, that is the lib - ration response in an ellipsoid, has also been addressed 
in the literature ( Chan et aLll2011 ; Zhang et «Lll20li ). When the ellipsoid is axisymmetric 
with respect to the rotation axis, our asymptotic analysis can be very easily adapted to 
find the zonal flow which will not be fundamentally different from the one in the sphere. 
By contrast, when the container is no longer axisymmetric along the rotation axis, the 
driving mechanism of the zonal flow may change as it may no longer be associated with 
nonlinear interaction of boundary layer solution but to nonlinear interaction in the bulk 
of pressure forced solutions. 

The two issues considered in this section are also discussed in the recent paper by 
Noir et all (|2012fh 



We would like t o than k M. A. Calkins for having provided the numerical data published 
in lCalkins et all (|2010l ) for the comparison done in section 5. This work has benefitted 
from discussions with M. Le Bars, D. Cebron, and F. Busse. We would also like to 
acknowledge discussions with N. Rambaux and C. Moutoux concerning the astrophysical 
applications of our results to Galilean satellites and exoplanet systems. 
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Appendix A. The nonlinear coefficients Ai, B\ , C\ and D[ 

The coefficients appearing in expressions (|4.11a .b) are given by 



A 1 = 



sin 



16 s_ 3 s + 3 



B, = - 



(1-i) 



(4s_ 2 (s_ + i s+) 2 s+ 2 cos 2 6 + (-s_ 6 + s_ 4 s+ 2 - 4i 

+ cos0(-4s_ 2 s+ 2 (.s_ 2 + s+ 2 ) 2 + (s_ 4 + s+ 4 ) sin 2 0)) 
7 (s_ 2 + s+ 2 )(s_ +is+)((s_ -is+) 2 -cos6»])sin 4 e», 



s_ 3 S+ 3 - S _ 2 S+ 4 



s+ 6 ) sin 2 i 



16s_ 2 s + : 
sin 2 cos , . . , 



16s+ : 



(4s+ cosfl + sin 2 9) , 



sin t)cosU , , . , 



4-5 = 



16s_ 2 

sin 2 6> 
16s_ 3 s + 2 
sin 2 ^ 



-As- cos9 + sin 2 



(s+ 2 + cos 6) (4cos0 (s_ - s + ).s+ 2 s_ 2 + sin 2 6»(s+ 3 + s- 3 )) 



3-(s_ 2 - cos 6») (4cos6»(s_ - s+)s+ 2 s_ 2 + sin 2 6>(s+ 3 + sJ)) 



I6s- 2 s 
B2 = -B3 = -B4 = -B5 = 0, 

Cl = (1 f ,~ i) 3 Sin 3 g ((s- 5 - S - 3 s + 2 - Ls_ 2 ,s+ 3 + i S+ 5 ) sin 2 I 



32s_ 3 s+ 



+4s_ (s_ - is+)(s_ 2 + s+)s+ cos 6) 



£>i 



16s -S -1 



° A = { \2s^f + e ( 4cos ^ s - - s +>+ 2s - 2 + sin2 d ( s + 3 + s - 3 )) > 
C 5 = - (1 ~ i) f na<51 (4cos0 - S+ ).s + 2 ,s_ 2 + sin 2 0(s + 3 + sJ)) , 

OZS+^S- 

C 2 = D 2 = C 3 = D 3 =D 4 = D 5 = 0. 



(Ala) 
(A 16) 

(Ale) 

(Aid) 

(Ale) 

(Alf) 
(A Iff) 

(Alh) 
(AH) 

(Alj) 

(A lk) 
(All) 



where 



s+ 



cos 6, 



(A 2) 
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Appendix B. The function T(p; ui) 

For bj > 2, we have the following expression for u) with p = sin9: 



J r (sm9; uj) 



(s + 3 + s_ 3 ) sin 4 6> (s_ - s+) cos 6> sin 2 6> 



32 s + 2 s. 



cost 



cos( 



s+ + cos( 



cos 9 



sm' 



9{s+ + V> 



cos ( 



+ 



32s+ 3 (cos 2 (9 + 4s + 4 , 

sin 2 9(s- + Vcos 9) 
32s_ 3 (cos 2 6» + 4s_ 4 ) 



- [2 s+ 2 - 2 s+ Vcos6> + cos 0] [sin 2 + 4 s+ 2 cos 0] 



[2 ; 



2 s_ Vcos9 + cos 0] [sin 2 9 - 4 s_ 2 cos 0] 



sm 



32 s_ 3 s+ 3 (s_ 2 + s+2) [s_ 2 + s+2 + 2 s + Vc~c~s0 + cos0] 2 
(sJ + S+ 2 )[s- 6 + 7 s+ 2 ,s_ 4 + 8 s_ 7 s+ 3 - 7 S _ 2 S+ 4 + 16 sJ s+ 5 

-s+ 6 + 4s+ 3 s_ 3 (-3 + 2 s+ 4 )] + [4s_ 6 .s+ - s_ 4 s+ 3 - 16s_ 2 s+ 5 - 3s+ 7 
+4.s_ 5 s+ 2 (-3 + 16s+ 4 ) + s+ 4 s_ 3 (-19 + 32s+ 4 ) + s_ 7 (-l + 32s+ 4 )] Vcos 6 



l (-ll + 4s+ 4 



+ [4s_ 8 s+ 2 - 3.s+ 6 + s_ 6 (l - 20.s+ 4 ) + s_ 

+4s+ 3 s_ 3 (-3 + 16s+ 4 ) + s_ 5 (-4s+ + 64s+°) + s_ 4 (s+ 2 - 20s+ e )] cos0 
-[4 s_ 7 s+ 2 + 4 s_ 6 s+ 3 + s+ 5 + 56 sJ s+ 5 + s_ 3 s+ 2 (3 - 76 s+ 4 ) 
+s_ 5 (1-8 s+ 4 ) + 3 s+ 3 s_ 2 (1-4 s+ 4 )] cos 3 / 2 9 

+ [-s_ 8 - 4s_ 6 s+ 2 - 4s+ 3 s_ 5 - 48s_ 4 s+ 4 + 60s + 5 s_ 3 + 20s_ 2 s+ 6 + s + 8 ] cos 2 



-[sJ - 4s_ 6 s+ + 8s_ 5 s+ 2 - lls_ 4 s H 



31s, 



20s_ 2 S+ 5 + 3 S+ 7 ] cos 5 / 2 



4s_ 



2 4 
S + S_ 



12s + 3 s_ 3 + lls_ 2 s + 4 + 3s+ 6 ] cos 3 



+ (s_ 5 + 3 s_ 3 s+ 2 + 3 s+ 3 s_ 2 + s+ 5 ) cos 7/2 6> 



+ - 



sm 



(Bl) 



where 



— + cos 9 and s_ 



--cos( 



(B2) 
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For u < 2, the above expression applies when p > p c = yX — w 2 /4. When p < p c , we 
have to use the following expression with p = sin 9 



J : (sm9; ui) = 



where 



sin 



32 



(Z+ + Vcosd) (4 1\ cos 9 + sin 2 9) {I- + V~cos9) (4 P_ cos 9 + sin 2 9) 
1% (2 1\ + 21+ Vcosd + cos 9) l 3 _ (2 1 2 _ + 2Z_ Vcostf + cos 9) 



4(l 3 _+l 2 _Vcos9+(l-+l+) cos 9 lism 2 9(l- + Vms9+(l++l-)/li 



l+l- (l- + V~c 



cos 91+ E (Z_ + Vcosl 



4^+4v / ^+((+ + L) cos9 l + sm 2 9(l + + V^9 + (l 3 _+ll)/l + cos9) 



Z_ Z + (Z + + vcosi 
1 



cos 6 1 1+ {1+ + \/cosl 



cos 6> ZrL /it (/- + 1+) (I- +1+ + Vcosd) 2 
^4Z 2 Z 2 (i_ - Z+) 2 cos 5/2 6> [1 + 2 (I_ + Z+) VW] 

+ cos 9 [I- 1+ {I- + 1+) (2 (2Z 2 - 31-1+ + 2l 2 +) sin 2 9 + Al 5 _l+ - 81^1 + Al-f+) 
+2 (/i + l\) sin 2 6»] + (Z_ + l+f {it + I*) sin 2 9 [(I- +l+) + 2 Vc^s~9] 



+ cos 3 / 2 9 [8 Z 2 Zi (Z 2 - Z 2 ,) 2 + {it + l 4 + + Z_Z+ (Z_ - l+f) s 



sm 



(B3) 



Z + = < /cos 



and Z_ 



cos9 



(B4) 



Appendix C. Stewartson layers 

In this section, the solutions in the different layers around p = a are provided. The 
analysis mainly follows the derivation given by IStewartson (1966). In all the layers, the 



functions ?p2 and X2 are related at leading order in E by a same set of equations which 
can be written in cylindrical coordinates as 



_ 2 * =E #» (cli) 

oz ap z 

In the two outer layers (see figure 0]), the function X2 remains at leading order in E 
independent of z, such that a simple relation between ip2 and X2 can be obtained by 
integrating (jC lb ) with respect to z: 



^ 2 (p, z ) = -E-^ + B(p), (C2) 



z d 2 X 2 
2 dp 

where B(p) is a function to be determined in each outer layer 



External outer layer E 1 / 4 

In the external outer layer of thickness E 1 / 4 , the function ^2 must be antisymmetric 
with respect to the equator z = 0. This implies B{p) — in this layer. Another relation 
is obtained by considering the problem close to the boundary layer at z = yl — a 2 . In 
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this relation between tp2 and X2 is deduced by writing (|4.18p close to sin 6 = a 

X2 (p) = -2 (1 - a 2 ) 1 / 4 E- 1 ' 2 MP, Vl-a 2 ) + F{a; to). 



Evaluating (IC 21) at z = yl — a 2 then leads to the equation 

<9 2 X2 X2--7 r (a;w) 



which immediatly gives 

X2O) = J"(ffl;w) - Ai exp 



dp 2 E 1 / 2 (1 - a 2 ) 3 / 4 ' 

(p-a)^- 1 / 4 



(1 - a 2 ) 3 / 8 



(C3) 
(C4) 

(C5) 



At this level, the constant Ai is, a priori, unknown. We shall see below that it is given 
at leading order in E by 



Ai = SX = X + - X~ = T{a\ lj) 
In this layer, the axial vorticity is given by 



4cj 2 



1 5x2 



1/4 6X(a,u)) 



■ exp 



{p-a)E- 1 ' A 
(1-a 2 ) 3 / 8 



Z2 p dp " a (1-a 2 ) 3 / 8 
while the function ip2 is given by the relation (|C 2j) with -B(p) = 0: 

(p- a) 



(1 - a 2 ) 3 / 8 ^ 



(C6) 



(C7) 



(C8) 



Internal outer layer E 2 ^ 7 
The solution in the Internal outer layer is obtained by considering the problem close to 
the two boundaries (at z = and z = \/l — a 2 ). Close to the boundary on the inner 
sphere (z = 0), we obtain a relation between ^2 and X2 by using (|4.20|) : 



Xa(p) = 2 [1 - (p/a) 2 ] 1/4 i?- 1 / 2 ^(p, 0) + 
Considering (|C 2j) at z = 0, we then obtain B(p) in this layer: 



4cj 2 



(C9) 



(CIO) 



If we now use the relation (|4.18p , we obtain from (IC 21) , applied close to z = yl — a 2 
with (IClOl) 



9 2 X2 X2-^(a) X2-a 2 a 2 /(4^ 2 ) a 1 / 4 



dp 2 (l-a 2 ) 3 / 4 \/£ [2(a-p)]V 4 v^yT^ 2 
This equation leads to the following solution 

x 2 (p) = ^ + £ 1/2 %f(d + 0(£ 1/M ), 

4 a; 2 

with 

a 



= 0. 



p = (a - p) 



2£ 2 (l-a 2 ) 2 



1/7 



(Cll) 



(C12) 



(C13) 



Libration-induced mean flow in a spherical shell 27 

and 

/2\ 1/7 (I - a 2 )- 5 / 56 
A 2 = -(-j L_^_« (a|W)| (014) 

where the function f(s) satisfies 

0-^T = O, f(0) = l, f(oo) = 0. (C15) 
An explicit expression for f (s) is provided by 

2(4/7) 11 /V/ 2 K 4/7 (8/7 S 7 / 8 ) 

f(s) = f(ii77) ' (c 16) 

where is a modified Bes sel function of order v and r(^) is the Gamma function 
( Abramowitz fc Stegunlfl965l). Note that f'(0) = -(4/7) 8 / 7 r(3/7)/r(ll/7) » -1.2246. 

The expressions (jC 6j) and (|C 14[) for Ai and A2 are such that X2 and d P X2, obtained 
from (IC 3P and (|CM)ll for p > a and p < a respectively, are continuous at p = a. As shown 
by IStewartsonl (|1966T ). this continuity condition comes from the condition of matching 
with the solution in the inner layer. It guarantees the continuity of the axial vorticity 
which is given in the Internal outer layer by 

^ = E ~ 1M Wfiw^m e ^ (C17 > 

Using (IC 21) . (jC 101) . (IC 1 1|) and (jC 12|) . an expression for tfj 2 can be obtained: 

^-^('-tct) ^.-^-^' ^ (C18) 

Inner layer E 1 ^ 3 

In the inner layer, %2 is constant and equal to a 2 a 2 /(4o; 2 ) . The function ^2 has by 
contrast a more complex structure. As shown by IStewartson ( 19661 ). ^2 = 0{E 19 / i2 ) in 
the inner layer. More precisely, it can be written as 

^ 2 ~£ 19/42 A 3 (a, W )I(i?,Z), (C19) 

with 

*= ( i- ( y«U - (c2o<,) 



and I(i?, Z) which satisfies 



<9 6 I a 2 i 

4^ = 0, (C21) 



with 

(-2i?) 1/4 1 - 1 when Z ->• for # < 0, (C 22a) 

(-2i?) 1/4 I- (1 -Z) when i? -)• -00, (C22&) 

I(R,Z = 1)=0, (C22c) 

Z = 0) = for R > 0, (C22d) 

I^0asi?^+oo. (C22e) 
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An integral form of I(R, Z) can easily be obtained from (|C 21j) using Fourier transform 
techniques as 

v ; 2V4 7T 7 sinh(fc 3 /2) k 3 / 4 ' 

The dependence of ip2 with respect to ui is only in the function A3 which is found to 

be 

a 3/28 

A 3 (a,u>) = - 26/7f , (Q) (1 _ q2)11/84 ■ (C24) 

The axial how u Z2 which can be deduced in the inner layer from tp2 by 



_ 1 di^ 2 

l z 2 



and 



T(3/4) f 00 sinh [fc 3 (l - Z)/2] 
J <* Z) = I sinh( fc 3/2) J C ° S ^ 8 " **> fcV4 dfc - 



(C25) 



a dp 
is then given by 

u Z2 ~ £ 5/12 A 4 (a,w) J{R,Z), (C26) 

with 

A 3 (a, a;) ^(0,0;) 

A4^,W; a(1 _ a 2)l/6 2 6/7 f /( 0)a 25/28( 1 _ a 2 ) 25/84' 



(C28) 



For the comparison with the numerical results, both leading-order and second-order 
approximations are used. Second order approximations can be obtained by considering 
the OiE 1 / 28 ) correction to X~ deduced from (|C1~2")) : 

X- = ^T+£ 1/28 A 2 f(0). (C29) 
This provides a new expression for 5X: 

5X - i + C£ ;i/28 ' ( C3 °) 

with 

(2V' 7 (l-a 2 )-V 56 f(0) 

C ~"W ^ ■ (C31) 

If we use (jC30)) in the expressions (jC6)) , ([UT4]) . (|C24| and (|C27| for the coefficients 
Ax, A 2 , A 3 , A4, we obtain approximations valid up to 0(E 1 ^ 14 ) corrections. It is these 
second-order approximations which have been used in figures IWb). UOf b) and lllf b). 



